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Abstract 

Climate  warming  is  expected  to  degrade  permafrost  in  many  regions  of  the 
world,  including  Alaska.  Degradation  of  permafrost  has  the  potential  to 
dramatically  affect  soil  thermal,  hydrological,  and  vegetation  regimes. 
Projections  of  long-term  effects  of  climate  warming  on  high  latitude 
ecosystems  require  a  coupled  representation  of  soil  thermal  state  and 
hydrological  dynamics.  Such  a  framework  was  developed  to  explicitly 
simulate  the  soil  moisture  effects  of  soil  thermal  conductivity  and  heat 
capacity  and  its  effects  on  hydrological  response.  The  model  is  the  result  of 
coupling  the  Gridded  Surface  Subsurface  Hydrologic  Analysis  (GSSHA) 
model  with  the  Geophysical  Institute  Permafrost  Laboratory  (GIPL)  model. 
The  GIPL  model  simulates  soil  temperature  dynamics,  the  depth  of  seasonal 
freezing  and  thawing,  and  the  permafrost  location  by  numerically  solving  a 
one-dimensional  nonlinear  heat  equation  with  phase  change.  The  GSSHA 
model  is  a  spatially  explicit  hydrological  model  that  simulates  two 
dimensional  groundwater  flow  and  one-dimensional  vadose  zone  flow. 
These  two  models  were  combined  by  incorporating  the  GIPL  model  into  the 
GSSHA  model.  The  GIPL  model  is  used  to  compute  a  soil  temperature 
profile  in  every  two-dimensional  GSSHA  grid.  GSSHA  uses  this  information 
to  adjust  hydraulic  conductivities  for  both  the  vertical  unsaturated  soil  flow 
and  lateral  saturated  groundwater  flow.  Test  case  results  indicate  that 
freezing  temperatures  reduces  soil  storage  capacity  thereby  producing 
higher  peak  discharges  and  lower  base  flow. 
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Unit  Conversion  Factors 
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By 

To  Obtain 
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0.02831685 

cubic  meters 

cubic  inches 

1.6387064  E-05 

cubic  meters 

cubic  yards 
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meters 

gallons  (US  liquid) 
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meters 

miles  per  hour 

0.44704 
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millimeters 

ounces  (US  fluid) 
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cubic  meters 

pints  (US  liquid) 
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cubic  meters 

quarts  (US  liquid) 

9.463529  E-04 

cubic  meters 

square  feet 
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square  meters 

square  inches 
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square  meters 
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square  meters 

square  yards 

0.8361274 

square  meters 

Yards 

0.9144 

meters 

1  Introduction 


Future  climate  change  scenarios  indicate  significant  increases  in 
temperature  are  projected,  especially  near  Earth’s  poles.  This  has  the 
potential  to  significantly  impact  regions  of  the  world  where  permafrost 
currently  exist,  near  the  poles  where  temperatures  are  projected  to  increase 
the  most.  Changes  to  permafrost  and  to  the  seasonally  frozen  soil  regime 
have  the  potential  to  significantly  alter  the  local  hydrology  and  resulting 
ecosystem.  More  generally,  the  soil-freezing  characteristic,  a  relationship 
between  unfrozen  water  content  and  temperature,  is  relevant  for  any  mass 
transfer  processes  in  frozen  porous  media.  To  better  understand  the  long 
term  effect  of  future  climate  scenarios,  especially  at  the  higher  latitudes, 
interaction  of  soil  thermal  state  and  hydrological  dynamics  is  significant. 
Thus,  a  couple  framework  was  developed  to  simulate  interactive  effects  of 
soil  thermal  and  hydrological  dynamics.  Gridded  Surface  Subsurface 
Hydrologic  Analysis  (GSSHA)  model  was  chosen  as  the  parent  code  in  this 
modeling  framework.  GSSHA  (Downer  et  al.  2006)  is  a  fully  distributed, 
physics  based  model  that  includes  the  ability  to  simulate  overland  flow, 
infiltration,  saturated  groundwater  flow,  evapotranspiration  (ET),  snow 
accumulation  and  melting,  as  well  as  many  other  physical  processes.  Past 
coupling  efforts,  for  example  coupling  of  subsurface  storm  drainage  and  tile 
drain  in  GSSHA  (Ogden  et  al.  2011;  Pradhan  et  al.  2009)  has  demonstrated 
GSSHA  ability  to  simulate  important  surface,  subsurface  runoff  generation 
processes  and  to  explicitly  represent  fully  coupled  hydrodynamics.  In  this 
present  framework,  the  Geophysical  Institute  Permafrost  Laboratory  (GIPL) 
model  (Jafarov  et  al.  2012;  Marchenko  et  al.  2008)  is  incorporated  into  the 
GSSHA  parent  model.  The  GIPL  model  simulates  soil  temperature 
dynamics  and  the  depth  of  seasonal  freezing  and  thawing  by  numerically 
solving  a  lD  quasi-linear  heat  equation  with  phase  change. 
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2  Purpose 

The  purpose  of  this  report  is  to  describe  the  development  of  the  coupled 
framework  to  explicitly  model  the  soil  moisture  effects  of  soil  thermal 
conductivity  and  heat  capacity  and  the  resulting  effects  on  hydrological 
dynamics.  The  report  consists  of  the  technical,  theoretical,  and  conceptual 
details  in  coupling  the  GIPL  permafrost  model  to  the  hydrologic  model 
GSSHA.  The  effects  of  seasonal  freezing  and  thawing  on  hydrological 
dynamics  are  demonstrated  by  applying  the  coupled  system  to  simplified 
test  cases. 

This  document  is  an  addendum  to  the  original  GSSHA’s  User  Manual 
(Downer  and  Ogden  2006)  as  it  describes  the  details  for  developments 
after  the  user’s  manual  was  completed.  Additional  information  on  GSSHA 
can  be  found  in  the  GSSHA’s  User’s  manual,  and  the  GSSHA  wiki 
(http://www.gsshawiki.com/Gridded  Surface  Subsurface  Hydrologic  Analysis!. 
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3  GIPL  Coupling  in  GSSHA 

The  basis  of  GSSHA  is  a  two-dimensional  (2D)  uniform  grid  used  for  both 
surface  and  subsurface  computations.  Point  computations,  infiltration,  ET, 
etc.,  are  performed  within  a  grid  cell  and  the  point  responses  are  inte¬ 
grated  to  get  the  system  response,  overland  flow,  lateral  groundwater  flow. 

GIPL  is  an  implicit,  finite-difference,  numerical  model  which  solves  the 
one  dimensional  (lD)  non-linear  heat  equation  with  phase  change.  The 
process  of  soil  freezing/thawing  is  treated  in  accordance  with  relationships 
between  the  soil  unfrozen  water  content  and  temperature.  A  special 
enthalpy  formulation  of  the  energy  conservation  law  makes  it  possible  to 
use  a  relatively  coarse  vertical  resolution  without  loss  of  latent  heat  effects 
in  the  phase  transition  zone.  The  mathematical  description  section  gives 
more  detailed  information  on  the  enthalpy  method.  In  the  2D  grid,  the  soil 
thermal  state  provided  by  GIPL  is  a  point  process,  and  is  solved  for  each 
cell  in  the  2D  grid. 

The  spatial  variability  of  land-surface  and  hydrodynamic  parameters, 
including  subsurface  soil  moisture  state,  are  included  in  the  GSSHA  model, 
and  made  available  to  GIPL  during  simulation,  Figure  1.  GIPL  uses  these 
values  to  update  the  thermal  state  of  the  soil  and  passes  this  back  to 
GSSHA.  GSSHA  uses  the  thermal  state  of  the  soil  to  determine  whether  the 
soils  are  frozen  or  unfrozen.  This  information  is  used  to  adjust  saturation 
levels,  hydraulic  conductivities,  and  saturated  groundwater  media  thickness 
used  in  water  flow  computations.  These  computations  produce  updated 
values  of  groundwater  level  and  soil  saturation  that  are  then  used  in  the 
GIPL  model  to  produce  new  thermal  state  profiles  in  each  grid  cell.  This 
change  of  information  continues  for  the  duration  of  the  simulation,  as 
depicted  in  Figure  2. 

3.1  GIPL  Mathematical  Model 

The  GIPL  numerical  model  solves  the  Stefan  problem  (Alexiades  and 
Solomon  1993,  Verdi  1994)  of  phase  change  which  is  the  problem  of 
thawing  or  freezing  via  conduction  of  heat.  The  enthalpy  formulation  is  used 
in  the  solution  of  Stefan  problem  in  GIPL.  The  core  of  the  GIPL  numerical 
model  is  based  on  the  lD,  vertical,  quasi-linear  heat  conductive  equation 
(Sergueev  et  al.  2003): 


ERDC  TR-13-15 


4 


Figure  1 GIPL  as  a  permafrost  component  in  GSSHA 
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Figure  2  Schematic  of  GSSHA  GIPL  coupling/linkage. 
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*  Computational  nodal  information  exchange  and  linkage  of  coupled  GSSHA  and  GIPL. 
The  linkage  include  (a)  GSSHA  soil  moisture  and  GIPL  temperature  linkage  (b)  GIPL 
temperature  and  GSSHA  hydraulic  conductivity  and  transmissivity  linkage. 
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qt  =v-(k(x,t)vt(x,  t))  (1) 

where  x  is  a  vertical  spatial  variable  which  ranges  between  xu,  upper  depth 
of  the  computational  unit,  and  xl,  lower  depth  of  the  computational  unit,  f 
is  temperature  and  r  is  time.  The  term  k{x,t)  is  a  thermal  conductivity 
(Wm^K"1);  H(x,t )  is  an  enthalpy  function. 

t 

H(x,t)  =  Jc(x,s)ds  +  L6(x,t )  (2) 


where  C(x,s)  is  volumetric  heat  capacity  (Mdim'sKr1),  6(x,t)  is  volumetric 
unfrozen  water  content  (%)  and  L  is  the  volumetric  latent  heat  of  freeze/ 
thaw  (MJm~3).  Equation  l  requires  boundary  and  initial  conditions.  The 
upper  part  of  the  domain  corresponds  to  the  air  layer  which  is  at  two  meters 
height  above  the  land  surface.  The  fictitious  domain  formulation  (Marchuk 
et  al.  1986)  allows  embedding  seasonal  snow  layer  into  the  current  air  layer. 
The  Dirichlet-type  boundary  condition  is  used  as  an  upper  boundary 
condition 


t{Xu>T)  =  tai r  (3) 

where  fair  is  a  daily  averaged  air  temperature.  The  geothermal  gradient  is 
set  at  the  lower  boundary: 


dt(xe,r) 

dx 


9 


(4) 


where  g  is  geothermal  gradient,  a  small  constant  (Km-1).  For  the  initial 
temperature  distribution,  an  appropriate  ground  temperature  profile 
based  on  the  point  location  is  used. 


t(x,T0)  =  f0(x) 


(5) 


The  formula  for  unfrozen  water  content  9{x,t)  is  based  on  empirical 
experiments  and  has  the  following  form: 


0(x,t) 


9M 


1, 

ale  —  t 


f  >  t, 
t  <  t, 


(6) 
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Parameters  a  and  b  are  dimensionless  positive  constants  (Lovell  1957),  c  is 
freezing  temperature,  while  r]  (x)  is  a  volumetric  soil  moisture  content.  The 
constant  t*  is  a  freezing  point  depression,  the  temperature  at  which  ice 
begins  to  form  in  the  soil,  so  that  there  is  no  ice  if  t  >  t* .  The  unfrozen  water 
content  9(x,  t)  varies  with  depth,  and  time,  and  is  dependent  on  soil  type 
and  hydrologic  forcing.  The  discritized  form  of  Equation  1  can  be  found  in 
Sergueev  et  al.  (2003)  and  Marchenko  et  al.  (2008).  A  detailed 
mathematical  description  of  the  model  and  numerical  solution  methods  can 
be  found  in  Nicolsky  et  al.  (2007). 

Required  input  data  include  climate  data,  snow  cover,  soil  thermal 
properties,  lithological  data,  and  vegetative  cover.  The  main  purpose  of  the 
model  is  to  generate  a  spatial  and  temporal  dataset  of  permafrost  distribu¬ 
tion  and  ground  temperature  dynamics  as  well  as  the  active  layer  thickness 
which  are  useful  in  a  wide  range  of  hydrologic,  ecological,  climatologic,  and 
socioeconomic  assessments  in  cold  regions. 

3.2  GSSHA  Subsurface  Processes 

3.2.1  Unsaturated  Zone  Model 

In  the  GSSHA  formulation  as  linked  to  GIPL,  a  lD,  vertical,  unsaturated 
model  rest  atop  a  2D,  lateral,  saturated  groundwater  flow  model,  as 
described  below.  GIPL  is  linked  to  GSSHA  in  both  of  these  domains. 

While  various  representations  of  the  unsaturated  zone  are  available  in 
GSSHA,  GIPL  is  linked  in  GSSHA  in  the  unsaturated  zone  with  the 
Richards’  Equation  (Richards  1931).  The  Richards’  Equation  is  a  general 
solution  of  saturated/unsaturated  water  movement  and  soil  moisture  and  as 
implemented  in  GSSHA  can  be  used  to  calculate  runoff  resulting  from  a 
variety  of  conditions,  including  infiltration  excess  and  saturation  excess 
mechanisms,  which  can  occur  simultaneously  in  different  areas  of  a 
watershed.  For  Richards’  Equation  there  is  no  requirement  that  the  runoff 
production  mechanism  be  known  a  priori  or  limited  to  one  type.  The 
GSSHA  model  uses  a  one-dimensional  finite-difference  solution  of 
Richards’  Equation  to  simulate  the  unsaturated  zone.  In  GSSHA,  the 
unsaturated  zone  is  linked  to  a  two-dimensional  finite-difference  represen¬ 
tation  of  saturated  groundwater  flow  (Downer  2002;  Downer  and  Ogden 
2004).  The  groundwater  solution  is  fully  coupled  to  surface  flows  using  a  lD 
implicit  finite  difference  solution  of  equation.  The  vadose  zone  controls  the 
flux  of  water  between  the  land  surface  and  groundwater  and  partitions 
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rainfall  into  runoff,  infiltration,  groundwater  recharge  and  ET.  In  GSSHA 
the  unsaturated  zone  below  each  overland  flow  cell  is  simulated  using  the 
one-dimensional  (vertical  direction)  head-based  form  of  Richards’  Equation 


c,„M 


dip 

d 

dr 

dz 

KJv) 


dip 


dz 


w  =  o 


(7) 


where  Cm  is  the  specific  moisture  capacity,  ip  is  the  soil  capillary  head 
(cm),  z  is  the  vertical  coordinate  (cm)  (downward  positive),  r  is  time  (h), 
KSOii(ip )  (cm)  is  the  effective  hydraulic  conductivity  and  W  is  a  flux  term 
added  for  sources  and  sinks  (cm  h1),  such  as  ET  and  infiltration.  The 
head-based  form  is  valid  in  both  saturated  and  unsaturated  conditions 
(Haverkamp  et  al.  1977). 

In  GSSHA  the  soil  column  is  subdivided  into  discrete  cells  and  Richards’ 
Equation  is  solved  using  a  cell-centered  implicit  finite-difference 
numerical  algorithm.  The  solution  scheme  is  central-difference  in  space 
and  forward  difference  in  time  and  is  thus  second-order  accurate  in  space, 
first-order  accurate  in  time.  Flux  updating  is  performed  to  ensure  mass 
balance  for  the  head  based  formulation. 


Variables  Ksoii  and  Cm  from  Equation  (7)  are  non-linear  on  the  water 
content  of  each  cell.  Unless  field  data  are  available,  the  Brooks  and  Corey 
(1964)  equations,  as  extended  by  Huston  and  Cass  (1987),  are  used  to 
calculate  Ksoii  and  Cm  based  on  the  water  content  of  the  cell.  One  exception 
occurs  when  there  are  saturated  cells  in  the  soil  column. 


3.2.2  Saturated  Groundwater  Model 

In  GSSHA,  the  2D  lateral  free  surface  water  flow  equation,  Equation  8, 
describes  the  movement  of  water  in  the  saturated  groundwater  zone.  The 
controlling  equation,  as  developed  by  Pinder  and  Bredehoeft  (1968): 


d_ 

dx 


dh) 

T  — 

XX 

OX 


d 

H - 

dx 


T  — 

v dy ) 


d 

H - 

dy 


_  dh 
T  — 

yx  dx 


d 

H - 

dy 


T  — 

mdy 


=  S^-  +  W{x,y,T)  (8) 
dr 


where  T is  the  transmissivity  (m2  s_1),  h  is  the  hydraulic  head  (m),  S  is  the 
storage  term  (dimensionless),  and  W is  the  flux  term  for  sources  and  sinks 
(m  s1).  It  is  assumed  that  off  diagonal  terms  are  not  important,  and  that 
transmissivity  can  be  expressed  as  the  product  of  the  hydraulic  conductivity 
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of  the  media  ( KSOii)  and  the  depth  of  the  saturated  media  (5).  Substituting 
surface  water  elevation  (Ews  =  h+ datum)  for  head,  the  free  surface  problem 
can  be  described  as  (Downer  2002) 

d  ,  OE  d  7  dE  9E 

—  Kb — —  H - Kb — —  =S — —  +  W(x,y,z)  (9) 

dx  '  dx  dy  dy  j  dr 

This  equation  is  solved  using  a  five  point  implicit  finite  difference  scheme. 
Solution  is  by  linear  successive  over  relaxation  (LSOR)  with  Picard 
iterations  on  both  T  (from  T  =  Kb)  and  S  (Downer  2002). 

3.3  Coupling  GIPL  to  GSSHA 

The  GIPL  model  is  a  standalone  permafrost  model  that  is  used  to  compute 
a  one-dimensional  (vertical)  soil  temperature  profile  over  time  using  static 
values  of  soil  moisture  at  daily  intervals.  As  implemented  in  GSSHA,  GIPL 
is  a  subroutine  that  is  used  to  compute  a  profile  of  soil  temperature  in 
every  2D  grid  cell,  including  time  varying  soil  moisture  and  groundwater 
levels  at  varying  time  intervals,  Figure  2.  To  accomplish  this  result  several 
tasks  were  performed: 

1.  The  original  Geophysical  Institute  Permafrost  Lab  (GIPL)  permafrost 
model  was  coded  in  FORTRAN.  This  FORTRAN  source  code  was 
converted  to  stand  alone  C/C++  source  code. 

2.  Originally,  GIPL  parameters  were  uni-dimensional  in  the  soil’s  vertical 
profile  but  are  lumped  in  the  horizontal  spatial  extent  of  application. 
Significant  effort  was  expended  to  make  all  the  GIPL  state  variables  and 
parameters  distributed  as  grid  based  or  permafrost  soil  type  based  before 
merging  the  C/C++  version  of  GIPL  into  GSSHA.  Thus,  the  uni¬ 
dimensional  limitations  of  GIPL  are  enhanced  into  multi-dimensional 
distributed  applicability  in  GSSHA  distributed  modeling  framework. 

3.  Originally,  the  GIPL  numerical  model  of  heat  transport  used  daily  or  larger 
time-steps.  As  implemented  in  GSSHA,  GIPL  can  have  any  time-step,  as 
specified  by  the  user.  The  default  time-step  is  the  infiltration  time-step, 
which  is  on  the  orders  of  seconds  or  minutes. 

4.  Several  thermo-hydrodynamic  formulations  and  modeling  concepts  are 
implemented  to  link  and  exchange  the  information  in  GIPL  and  GSSHA, 
as  described  in  Section  3.3.1 
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3.3.1  Linking  GIPL  and  GSSHA  Computational  Nodes 

In  GSSHA,  the  unsaturated  zone  is  subdivided  into  computation  nodes 
that  cover  the  time  varying  unsaturated  zone,  which  rest  on  top  of  the 
saturated  groundwater.  The  unsaturated  zone  is  divided  into  four  regions, 
corresponding  to  the  A,  B,  and  C  soil  horizons,  as  well  as  the  groundwater 
media.  The  saturated  groundwater  computation  in  2D  lateral,  so  there  is 
no  distribution  in  the  vertical  direction  in  the  saturated  zone.  The  location 
of  the  water  table  may  be  anywhere  from  or  above,  the  land  surface,  to 
hundreds  of  meters  below  the  soil  surface.  In  GIPL,  the  soil  column 
extends  down  to  where  the  lower  boundary  condition  is  considered  valid, 
Equation  (4),  very  deep  within  the  permafrost  region.  This  may  extend 
1000  or  more  meters  below  the  land  surface.  Because  of  the  differences  in 
domains  and  requirements  for  solution,  in  the  coupled  framework,  the 
GIPL  computational  nodes  are  independent  in  terms  of  location  of  the 
GSSHA  infiltration  scheme  computational  nodes.  Thus,  a  user  of  the 
permafrost  model  in  GSSHA  does  not  need  to  spend  time  matching  the 
vertical  discritization  in  GIPL  and  GSSHA  infiltration  numerical  schemes. 
The  linkage  of  computational  nodal  discretized  information  from  GIPL  to 
GSSHA  and  vice-versa  is  shown  in  Figure  2. 

3.3.2  Linking  GIPL  Soil  Thermodynamics  with  GSSHA  Soil  Moisture 
Hydrodynamics 

The  main  information  being  passed  during  a  GSSHA/GIPL  coupled 
simulation  is  that  GSSHA  provides  updated  soil  moisture  profiles  to  GIPL, 
which  is  used  to  adjust  the  thermal  capacity  of  the  soil,  and  that  GIPL 
provides  updated  soil  temperatures  to  GSSHA,  which  is  used  to  adjust  the 
hydraulic  conductivity  of  the  soil,  both  vertically  and  laterally.  The 
infiltration  time  step  with  updated  soil  moisture  of  GSSHA  is  used  to 
update  the  soil  temperature  profile  during  a  GIPL  time  step.  The  linkage  of 
GSSHA  soil  moisture  update  to  GIPL  thermal  update  is  shown  in  Figure  2. 

3.3.3  Linking  GIPL  Soil  Temperature  and  GSSHA  Hydraulic  Conductivity 

In  the  unsaturated  zone,  the  temperature  effects  are  accounted  for  by 
adjusting  the  value  of  relative  saturation  (Se)  and  using  that  adjusted  value 
to  adjust  the  vertical  hydraulic  conductivity  of  the  soil,  based  on 
temperature  and  saturation. 
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3.3.3. 1  Estimation  of  the  Relative  Saturation: 

The  relative  fraction  of  liquid  water  of  the  total  soil  moisture,  Se  is  defined 
as  follows  (Schulla  2012): 


SE  = 


l  +  (oc|l.22f|)" 


fort  <0°C 


(10) 


where  n,  m,  and  a  are  the  Van-Genuchten-Parameters  as  used  in  the 
Richards  Equation;  t  is  soil  temperature  in  °C.  For  temperatures  above  o°C, 
Se  is  always  1;  whereas  for  temperatures  below  -io°C  the  value  of  Se  is 
assumed  to  be  o.  The  latter  assumption  is  to  reduce  computational  burden. 
Fine  soils  may  still  contain  liquid  water  below  -io°C,  but  the  fraction  is 
negligible  to  any  thermal  and  hydraulic  transport  process  at  the  timescales 
of  model  applications  (Schulla  2012). 

3. 3. 3. 2  Linking  GIPL  Soil  Profile  Temperatures  and  GSSHA  Effective 
Hydraulic  Conductivity: 

The  effect  of  reduced  effective  saturation  in  the  soil  due  to  freezing  is 
reduced  hydraulic  conductivity.  An  effective  hydraulic  conductivity  is 
computed  for  the  entire  soil  matrix,  including  both  the  unfrozen  and 
frozen  fractions. 

In  the  unfrozen  portion  of  the  soil  an  exponential  response  in  effective 
hydraulic  has  been  measured  for  freezing/thawing  mineral  and  organic  soils 
(Zhang  et  al.  2010).  Accordingly,  the  exponential  function  is  applied  to 
calculate  the  effective  hydraulic  conductivity  where  the  hydraulic  conduc¬ 
tivity  K  at  a  given  temperature  (t)  is  a  function  of  hydraulic  conductivity  of 
the  unfrozen  soil  and  the  effective  saturation  Se  is  as  follows: 

(t )  =  es-ln  (K,  (©))■ +  (1  ■ - . SE  )Kf  (11) 


where  Ksoii(t)  is  the  effective  hydraulic  conductivity  in  m/s;  Kt  is  the 
hydraulic  conductivity  for  Se  =  1  and  Kf  is  the  frozen  hydraulic  conductivity 
( Se  =  o).  In  practice,  the  contribution  from  the  frozen  portion  of  the  soil, 

( i-Se )  Kf  is  quite  small  and  is  often  neglected.  The  linkage  of  GIPL  thermal 
nodal  information  and  GSSHA  hydraulic  conductivity  nodal  information  is 
shown  in  Figure  2. 
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3.3.4  Linking  Soil  Heat  Transfer  Effect  on  Effective  Groundwater 
Transmissivity 

As  described  in  Section  3.2.2,  the  transmissivity  (T)  of  the  saturated  media 
is  the  product  of  the  lateral  hydraulic  conductivity  ( K)  and  the  thickness  of 
the  saturated  media  (6).  As  GIPL  provides  temperatures  along  the  entire 
soil  profile,  including  both  saturated  and  unsaturated  (see  Section  3.2.1) 
zones  and  the  soil  profile  representations  in  GIPL  and  GSSHA  are  also 
linked  (see  Section  3.3.1)  it  is  possible  to  determine  what  portions  of  the 
saturated  zone  are  frozen,  and  which  are  not.  Since  the  representation  of 
the  saturated  groundwater  media  in  GSSHA  is  a  2D  lateral  free  surface 
equation  (see  Section  3.2.2)  the  frozen  layers  cannot  be  explicitly 
represented  as  such  in  the  GSSHA  model,  but  the  saturated  thickness  (the 
depth  of  water  flowing  in  an  unsaturated  flow  cell)  can  be  adjusted.  In  the 
coupled  framework,  the  thickness  of  the  effective  saturated  media,  b,  is 
computed  by  identifying  the  unfrozen  sections  of  the  soil  profile  and 
accumulating  those  unfrozen  layers.  This  effective  saturated  thickness  is 
used  to  compute  T  in  the  GSSHA  groundwater  sub-routine  while  updating 
the  groundwater  heads. 

The  depth  of  the  unfrozen  saturated  media  in  GSSHA  is  determined  by 
searching  the  GIPL  nodes  corresponding  to  the  saturated  media  depth.  If 
the  GIPL  node  is  within  the  saturated  media  and  is  not  frozen,  the  thickness 
of  that  node  is  added  to  the  saturated  media  depth.  The  search  for  frozen 
temperatures  in  the  GIPL  computational  nodes  is  made  proceeding  from 
the  top  of  the  soil  layer  to  the  deepest  computational  node  depth.  If  the 
nodal  soil  temperature  is  positive,  the  corresponding  computation 
block/ slice  dimension  is  added  to  the  saturated  media  depth.  Otherwise  the 
nodal  soil  temperature  is  negative,  and  the  corresponding  computation 
block/ slice  dimension  is  not  added  to  the  saturated  media  depth. 

The  applicable  block/ slice  dimension  is  added  to  the  saturated  media 
thickness  only  if  the  corresponding  GIPL  nodal  elevation  is  within  the 
limit  of  GSSHA  groundwater  table  elevation  and  the  GSSHA  aquifer 
bottom  elevation. 

In  this  top-to-bottom  approach  where  ]’  is  a  GIPL  node  number,  if  the 
node  above  or  below,  the  j-f  nodal  temperature  is  frozen  and  the  f  nodal 
temperature  is  unfrozen,  the  thickness  of  the  saturated  media  is  defined  as 
the  interpolated  unfrozen  depth  between  the  frozen  node  and  unfrozen 
node.  This  avoids  the  overestimation  of  the  effective  saturated  depth. 
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Once  the  effective  saturation  depth  is  calculated,  local/grid  based  GSSHA 
groundwater  transmissivity  is  defined  as  the  following: 


T 


_  k  * 

groundwater 


B 


effective 


(12) 
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4  Test  Cases 

The  test  case  example  in  this  section  illustrates  modeling  a  permafrost 
active  area  with  GIPL  coupled  in  GSSHA.  The  simplified  example  is 
conceptual  but  the  permafrost  parametric  values  represent  Alaskan 
woodland  and  tundra  ecosystem  sites  in  permafrost  active  regions.  This 
example  project  includes  surface  subsurface  runoff  where  infiltration  and 
groundwater  components  are  turned  on.  The  soil  moisture  and  soil 
physical  state  is  defined  by  the  Richards  Equation. 

4.1  Example  Model 

Figure  3  shows  the  test  case  example  model  having  10x10  building  blocks. 


Figure  3.  A  test  case  10x10  example  project  of  coupled  GSSHA  and  GIPL, 
where  the  permafrost  parametric  values  represent  woodland  and  tundra 
ecosystem  sites  in  permafrost  active  Alaskan  regions. 
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4.1.1  Permafrost  Boundary 

The  entire  test  case  is  regarded  as  a  permafrost  active  zone.  The 
permafrost  soil  properties  as  defined  in  Table  l  are  from  an  Alaskan 
woodland  and  tundra  ecosystem  permafrost  site. 


Table  1.  Permafrost  parametric  value. 


Description 

unit 

value 

Volumetric  soil  water  content 

Fraction  of  1 

0.87 

Volume  of  unfrozen  water 

Fraction  of  1 

0.11 

A-parameter  of  unfrozen  water 

- 

0.034 

B  parameter  of  unfrozen  water 

- 

-0.32 

C  parameter  of  unfrozen  water 

- 

0.0 

Soil  thermal  conductivity  thawed 

W  m  1  k 1 

0.0201 

Soil  thermal  conductivity  frozen 

W  m  1  k1 

0.0551 

Volumetric  heat  capacity 

JirAm-im-ik-1 

2800 

4.1.2  Initial  Condition 

Figure  4  shows  the  initial  temperature  condition  which  is  from  an  Alaskan 
woodland  and  tundra  permafrost  location. 


Figure  4.  Soil  temperature  profile  as  an  initial  condition  for  the 
thermodynamics  numerical  simulation. 
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4.2  Model  Results 

Figure  5  shows  simulated  soil  temperature  profile  extracted  from  the  time 
series. 


Figure  5.  Soil  temperature  profiles 


Figure  5  shows  the  vertical  soil  temperature  at  computational  nodal  points 
and  Figure  6  shows  the  corresponding  depths  for  the  nodal  points. 


Figure  6.  Depth  information  of  the  computational  nodal  number. 
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Permafrost  is  separated  from  the  atmosphere  by  a  boundary  layer 
consisting  of  the  active  layer.  The  active  layer  transmits  heat  to  and  from 
permafrost.  The  top  of  permafrost  is  at  the  base  of  the  active  layer.  An  active 
layer  is  the  soil  column  that  experiences  normal  freezing  and  thawing 
during  the  seasons.  Figure  5  shows  the  freeze  cycle  of  the  active  layer.  In 


ERDC  TR-13-15 


16 


Figure  5,  the  temperature  profile  in  the  active  layer  shows  decreasing  values 
while  moving  from  one  time  series  to  another.  This  decrease  in  the 
temperature  values  lead  to  decrease  in  free  water  content  in  the  frozen  soil 
thereby  decreasing  the  hydraulic  conductivity  of  the  soil.  Computational 
nodes  were  closer  in  the  active  layer  which  is  shown  in  Figure  6. 

Figure  7  shows  the  soil  temperature  at  various  depths.  Figure  7  shows  that 
the  air  temperature  has  the  most  significant  influence  in  the  near  surface 
soil  layer.  As  the  soil  layer  depth  increases,  air  temperature  influence  in 
soil  thermo-dynamics  is  diminished  along  with  the  increase  in  the  time  lag 
influence. 


Figure  7.  Time-series  of  temperature  at  various  depths. 


Time  (hrs) 


Figure  8  shows  the  change  in  effective  hydraulic  conductivity  due  to  frozen 
soil  condition.  The  effective  hydraulic  conductivity  decreases  with 
increasing  fraction  of  ice,  i.e.  a  decreasing  Se  value  in  Equation  (10)  and 
Equation  (11).  The  effective  hydraulic  conductivity  changes  with  several 
orders  of  magnitude  as  the  soil  freezes/thaws  which  is  defined  by  Equation 
(11)  with  an  exponential  response  in  effective  hydraulic  conductivity.  Figure 
8  shows  the  decrease  in  soil  effective  hydraulic  conductivity  from  1  cm  hr1 
to  almost  zero  when  the  active  layer  soil  column  started  to  freeze.  This 
simulation  result  of  the  change  in  effective  hydraulic  conductivity  due  to  soil 
freezing  condition  agrees  with  the  fact  that  the  effective  hydraulic  conduc¬ 
tivity  changes  by  several  orders  of  magnitude  as  the  soil  freezes/thaws. 
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Figure  8.  Hydraulic  conductivity  under  active  permafrost  soil  layer. 


Figure  9  shows  the  comparison  of  GSSHA  simulated  discharge  with  and 
without  the  permafrost  model. 

4.3  Discussion 

Stream  response  to  storms  differs  for  permafrost-free  and  permafrost- 
affected  slopes.  Stream  flow  from  watersheds  underlain  with  a  large 
proportion  of  permafrost  responds  rapidly  to  precipitation,  with  rapid 
rising  and  falling  limbs  of  storm  hydrographs  (Quinton  and  Carey  2008). 
The  simulation  result  with  permafrost,  shown  in  Figure  9,  had  rapid 
response  to  the  precipitation  event  with  increased  peaks  than  that  without 
the  permafrost.  This  verified  that  the  model  produced  desired  and 
expected  output. 

In  contrast,  for  streams  with  little  or  no  permafrost,  precipitation  can 
percolate  into  soil  layers  resulting  in  enhanced  connectivity  between  the 
surface  and  ground  water  storage  regimes  and  more  soil  pore  water 
storage.  Because  of  this  enhanced  connectivity  and  soil  pore  water  storage 
capacity,  Figure  9  shows  reduced  peak  discharge  and  runoff  volume  for  the 
permafrost-free  simulation  in  comparison  to  the  result  with  permafrost. 

All  those  simulation  results  presented  in  Section  4.2  show  that  GSSHA 
coupled  with  GIPL  could  serve  as  a  valuable  tool  for  long-term  simulation 
and  prediction  of  interactive  effects  of  frozen  soil  hydrological  dynamics  in 
permafrost  regions. 
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Figure  9.  Hydrograph  with  and  without  active  permafrost. 
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5  Summary 

A  coupled  framework  was  developed  for  simulating  the  interaction 
between  soil  temperature,  including  permafrost,  and  hydrology,  by 
incorporating  the  soil  temperature  and  permafrost  model  GIPL  into  the 
distributed,  physics  based  hydrologic  model  GSSHA.  The  report  describes 
the  numerical  considerations  in  linking  the  GIPL  thermo-dynamic  model 
into  GSSHA’s  hydrodynamic  modeling  framework.  GSSHA  hydrodynamics 
include  soil  moisture  saturation  feedback  in  the  vadose  zone  and  the 
corresponding  soil  ice  content  effects  on  hydraulic  conductivity  and 
transmissivity. 

The  coupled  model  was  demonstrated  on  a  contrived  watershed,  around  a 
previous  GIPL  test  site.  The  coupled  simulation  results  show  that  the  effect 
of  soil  thermal  properties  obtained  from  GIPL  play  a  significant  role  in  the 
GSSHA  hydrological  dynamics  and  vice  versa. 
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